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The energy level statistics of the Hubbard model for LxL square lattices 
(L = 3, 4, 5, 6) at low filling (four electrons) is studied numerically for a wide 
range of the coupling strength. All known symmetries of the model (space, 
spin and pseudospin symmetry) have been taken into account explicitly from 
the beginning of the calculation by projecting into symmetry invariant sub- 
spaces. The details of this group theoretical treatment are presented with 
special attention to the nongeneric case of L = 4, where a particular compli- 
cated space group appears. For all the lattices studied, a significant amount 
of levels within each symmetry invariant subspaces remains degenerated, but 
except for L = 4 the ground state is nondegenerate. We explain the remain- 
ing degeneracies, which occur only for very specific interaction independent 
states, and we disregard these states in the statistical spectral analysis. The 
intricate structure of the Hubbard spectra necessitates a careful unfolding pro- 
cedure, which is thoroughly discussed. Finally, we present our results for the 
level spacing distribution, the number variance T?, and the spectral rigidity 
A3, which essentially all are close to the corresponding statistics for random 
matrices of the Gaussian ensemble independent of the lattice size and the 
coupling strength. Even very small coupling strengths approaching the inte- 
grable zero coupling limit lead to the Gaussian ensemble statistics stressing 
the nonperturbative nature of the Hubbard model. 

PACS numbers: 75.10Jm, 74.20.-z, 05.45. -hb 



I. INTRODUCTION 

The behavior of strongly correlated electronic systems remains a central problem in 
contemporary condensed matter physics. Several years of intense studies have made it clear 
that the necessary theoretical skills and tools to deal with strongly correlated fermion systems 
are lacking (see e.g. the recent reviews by Dagatto and Lieb [Q). Many exotic schemes 
have been invented to accommodate a suitable theoretical framework, but the development 
of a predictive general theory does not seem to be in sight. In this state of affairs the 
importance of performing numerical calculations of the ground state properties and the 
energy spectrum of a given many-body Hamiltonian has grown. Computational results can 
lead to the acceptance or rejection of the proposed analytical models, and they can guide 
the development of new analytical approaches. 

In the Hubbard model and related models one important parameter is the filling v. Much 
work has been devoted to the high density case near half filling, since it is believed to be 
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relevant for high-temperature superconductivity |T|J^, but also the low filling regime is of 
interest, e.g. it plays an important role in theoretical studies of the breakdown of Fermi 
liquid theory in 2D [Q. In this paper we present a numerical study of the two-dimensional 
Hubbard model at low filling regime v < 0.25, a regime where the calculation is tractable. 
It is natural to choose four particles as a generic case close to the simple two-particle case. 
The coupling strength is used as a perturbation parameter, and we address the question of 
universality in the response of strongly correlated electron systems to this perturbation [Q. 

We describe an efficient method which allows for numerical calculations of the exact 
energy spectrum. This method can relatively easily be extended to calculations of various 
Green's functions and spectral functions well suited for the study of low lying excitations 
and the corresponding coherent part of the spectral densities, a topic we will study in 
forthcoming work. Here, we rather study the statistical properties of the typical high energy 
excitations which is related to the incoherent background of the typical spectral functions. 
More specifically, we study the statistical properties of the spectra within the framework 
of random matrix theory (RMT). RMT was developed for the study of neutron scattering 
resonances in nuclear physics in the fifties and sixties p, but it has since been applied to a 
wide range of problems in many areas of physics (e.g. studies of conduction fiuctuations 
microwave eigenmodes [^], and acoustical properties of solids [jTU[) and mathematics 
(e.g. studies of the distribution of the zeros of the Riemann zeta function [0). Moreover, 
RMT has also been applied to various types of matrix ensembles like Hamiltonian matrices 
1^ (as in our case), scattering matrices PJl^, as well as transfer matrices |13| and Glauber 



matrices of statistical mechanics models. Recently, RMT has been employed in the 
study of strongly correlated electronic systems. Examples are studies of the 2D tJ-model 
p!5| , 2D tight-binding models [IC], the ID Bethe chain [17|, and the ID Luttinger liquid 



TE| . The work presented here with emphasis on mathematical and numerical methods is 



an extension of this line of research. Preliminary results of our work have been published 
elsewhere JT^j . 

There are basically two ways of applying RMT. One way is to model a relevant matrix 
of the given physical system with a matrix drawn from a suitable random matrix ensemble 
and subsequently calculate average properties of the system by averaging over the random 
matrix ensemble according to RMT. The other way, which we employ here, consists simply 
of characterizing the spectrum of a given physical system by comparing various statistical 
properties of the spectrum with the corresponding properties calculated within one of the 
few universal statistical matrix ensembles of RMT. Which of these ensembles describing 
properly a physical situation depends on the symmetries of the system. The given spectrum 
which one analyzes is of course deterministic, but statistical properties are given to it by 
considering quantities like for example the distribution of the energy spacings where all 
but one of the spectral variables have been integrated out. This is like pseudo-random 
number generators which are perfectly deterministic and nevertheless have many properties 
in common with random sequences. 

The paper is organized as follows. In Sec. II we introduce the Hubbard model and 
the corresponding Hilbert space. In Sec. HI we introduce the RMT quantities used in 
the characterization of the spectra, and we discuss in detail the special spectral unfolding 
technique necessitated by the intricate nature of the Hubbard spectra. In Sec. IV the 
entire symmetry group consisting of space, spin and pseudospin symmetry of the model is 
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studied, and all the corresponding projection operators are calculated. In Sec. V the model 
is diagonalized numerically and we study the raw spectrum, in particular the ground state 
and some unexpected remaining degeneracies higher in the spectrum. In Sec. VI we present 
the result of the spectral statistics analysis of the model, and finally, in Sec. VII we discuss 
the results and conclude. Appendices A - D contain mathematical details. 



II. THE HUBBARD MODEL 



Throughout this paper we study the simple one-band LxL square lattice Hubbard model 
with periodic boundary conditions containing a nearest-neighbor hopping term -tT and an 
on-site interaction term UU: 

L2 1,2 



H 



-tT + UU = -t c]^Ci„ + ^iT^a 



(1) 



where c]^ and the creation operator and the number operator, respectively, for an 

electron on site i with spin a. No disorder is present in the model. Below half filling the 
dimension Nh of the Hilbert space grows rapidly as a function of L and the number Ne of 
electrons occupying the lattice, so we have confined ourselves to the low filling case of only 
four electrons, whereas we let L vary. Moreover, without loss of generality we always work 
in the sector where the z component 5*2 of the total spin is zero; the other Sz sectors can be 
reached by use of the spin ladder operators S+ and S^, which commute with the Hamiltonian. 
The 5*2=0 sector is the largest of the spin sectors, and it has Nh = [L'^{L'^ — l)/2]^ which 
for L = 3, 4, 5, and 6, the lattice sizes studied here, yields 1296, 14400, 90000 and 396900, 
respectively. The corresponding fillings u = Nf./2L'^ are 0.22, 0.13, 0.08, and 0.06. 



(a) 



(b) 





FIG. 1. Two four-electron states with 5*2=0 are shown: in (a) the zero-pair state |5,6;2, 14) 
of the 4x4 square lattice, and in (b) the one-pair state |5, 16; 13, 16) of the 5x5 square lattice. 
The convention of the ket-notation is explained in the text. The numbering of the lattice sites is 
the one employed in our computer calculations. The signs of the 4x4 lattice correspond to the 
bi-partition of the lattice, which is only possible for even site lattices. 



In the occupation number basis we label the states as follows ||2^ , where a and b are two 
lattice sites occupied with spin-up electrons and c and d with spin-down electrons: 
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|a, b; c, d) = cl^cl^cl^cl^\vac). (2) 

Two explicit examples of such states as well as the lattice site enumeration are shown in 
Fig. |l|. In our work we make sure that a < b and c < d, and we have ordered the basis states 
such that the state \xi, X2; x^, X4) comes before the state \x[,X2; x'^, x'^) if Xi < x[, where i is 
the first position encountered where Xi and x'^ are different. If during a calculation a state 
is encountered with a > b and/or c > d, the necessary exchange operations including sign 
changes are performed to restore it. 

Finally we note, that for even L a bi-partition of the lattice is possible. A given lattice 
site a can be identified by a set of cartesian coordinates a = (01,02) simply counting the 
position in the lattice (thus e.g. site = (0,0) and 9 = (1,2) in Fig. |l]a). Each site a can 
then be assigned with a sign 6a = (— l)"i+"2 _ exp(z7r-a), where tt = (7r,7r). 

III. RANDOM MATRIX THEORY 

Within random matrix theory (RMT) one can study several statistical ensembles of ma- 
trices. Three important examples are the diagonal ensembles, the Gaussian ensembles for 
Hermitian matrices as e.g. Hamiltonians, and the circular ensembles for unitary matrices as 
e.g. scattering matrices. Since a main object of this work is to characterize the spectrum 
of the Hubbard model within the framework of RMT, we are led to use the diagonal and 
the Gaussian ensembles of square matrices. The first ensemble is the ensemble of diagonal 
matrices D with statistically independent diagonal elements Da drawn from the same dis- 
tribution P{Dii). This ensemble describes situations where the eigenvalues are essentially 
independent, which empirically has been found to be the case for integrable models. One 
characteristic feature of such systems is a "soft" spectrum with large probabilities of having 
levels close together described by the Poisson (exponential) distribution. The three others 
ensembles - denoted the Gaussian orthogonal ensemble (GOE), Gaussian unitary ensemble 
(GUE), and Gaussian symplectic ensemble (GSE) - are defined by requirering statistical 
independence of the matrix elements of a matrix H and invariance of the probability distri- 
bution P{H) in matrix space under one of the three canonical similarity transformations, the 
orthogonal, the unitary, and the symplectic transformation respectively [§. Which ensemble 
to choose depends on the symmetry of the system. In fact, the ensembles are universal in 
the sense that no details of the physical system plays any role, only knowledge of the global 
symmetry is needed. The three Gaussian ensembles are found to describe situations of very 
complex or chaotic systems, and one characteristic feature of such systems is a "rigid" spec- 
trum with eigenvalue repulsion. In this work we need only to treat GOE, which is found for 
systems with preserved time-reversal symmetry. 

To perform a meaningful RMT analysis one has to sort the spectrum in symmetry sec- 
tors corresponding to the symmetry group of the Hamiltonian since the symmetry invariant 
subspaces are orthogonal to one another. Each such symmetry invariant subspace is charac- 
terized by a specific set of quantum numbers, and the RMT analysis is performed on sets of 
eigenlevels having the same quantum numbers. In Sec. |3 we treat the complete symmetry 
group of the Hubbard model and calculate the corresponding projection operators of the 
symmetry invariant subspaces. As illustrated in Sec. ^ significant errors are introduced in 
the analysis if some of the symmetries are neglected. 
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Another caveat in the RMT analysis of finite spectra is the notion of mixed phase space. 
As a function of some external parameter a given system can be driven from integrability (the 
diagonal ensemble) to chaos (one of the Gaussian ensembles) or from one type of symmetry 
(say GOE) to another (say GUE). If the spectrum of such a system is studied in the middle 
of the transition the spectrum is a mixture of two or more components each described by 



one of the random ensembles [^. In the thermodynamic limit usually only one component 
survives, but for finite spectra this mixing calls for a further sorting of the spectrum within 
each symmetry invariant subspace. For the Hubbard model using the coupling strength U /t 
as the external parameter such a situation does in fact arise as the analysis presented in 
Sec. |V B| shows. 

When the necessary sorting of the spectrum has been performed the RMT analysis can 
begin. The first step is to unfold the spectrum. 



A. Unfolding the spectrum 

Naturally, it is necessary to make some kind of transformation or normalization of the 
spectrum of any given physical system to be able to make comparisons with the universal 
and dimensionless results of RMT. This operation is called the unfolding. By local rescaling 
of the spectrum with the local average level spacing the unfolding transforms the actual 
energies Ei into dimensionless "unfolded energies" with a local density of one. Thus, by 
unfolding one subtracts the regular slowly varying part of the spectrum and consider only 
the fiuctuations. It amounts to compute from the actual integrated density of states N[E), 

N{E) = r -^(e - Ei) de = Yl 0{e - E,), (3) 

% I 

an averaged integrated density of state N{E). The unfolded energies Ei are then given by 

e^ = N{E,). (4) 

The notions of "local density" and "averaged density" are not mathematically rigorous. For 
some systems there exists natural unfolding procedures. For example, it is a rigorous result 
that the density of states for a N x N random matrix approaches a semi-circular form for 
N ^ oo, thus for any given finite random matrix one simply uses the limiting density of 
states as the average density. In billiard systems one can make a Laurent series expansion 
in a/E of N{E) and obtain the Weyl law for N{E) by truncating the series after a finite 
number of terms, each of which has a physical interpretation. In our case, no such natural 
choices exist for N{E). We have therefore used several methods to unfold the spectra. Each 
of these methods has a free parameter, but there is no unique prescription of how to choose 
it. The best criterion is the insensitivity of the final result to the method employed and to 
reasonable variations of the free parameter. 

The first method is polynomial interpolation. It can be a simple linear interpolation or 
running average where N{Ei) is found as a linear fit of N{E) in an interval containing r 
levels on each side of the level Ef, the free parameter is then the parameter r. It can also 
be higher order polynomial interpolation, e.g. in the form of interpolating between several 
linear interpolations - a method we used in our work on statistical mechanics models ]r5|,p2 . 
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The second method |^ defines Ei = Si^i + di/dn, where dk is the fc'th smallest spacing 
to Ei. Here n becomes the free parameter. This method works also for complex eigenvalues 
of non-Hermitian matrices. 

The third method is Fourier broadening of the step-functions 6{e — E^) in Eq. @. The 
Fourier transforms from the energy domain to the time domain of the step-functions are 
found. In the following back-transformation yielding N{E) only the slow time components 
are kept. Choosing a cut-off r beyond which all Fourier components are set to zero yields 
9{e - Ei) Si[(e + E* - Ei)T]/7i - Si[(e - E* + Ei)T]/Tr, where Si(x) is the sine integral and 
E* is an energy slightly larger than the largest energy in the spectrum to be unfolded. The 
free parameter is r, and by choosing 1/r to be of the order of the mean-level spacing a good 
N{E) is obtained. 

The fourth method is Gaussian broadening of the delta-functions 6{e — Ei) in Eq. 
leading to the following expression for N{E): 




FIG. 2. Unfolding of the spectrum using the Gauss broadening with variable width choosing 
a = 4. Shown is a part of the Hubbard spectrum of the invariant subspace {R, S) = (6, 0) around 
E/t ~ 1.1 for L = 5 and U/t = 1. The positions of the levels Ei are marked by short vertical lines. 
For four levels marked by long vertical lines we show the actual broadened Gaussians. Note how 
the widths of the Gaussians change as the local density of states changes, and note how the width 
of the Gaussian centered near the gap ignores the states beyond the gap. Also shown are the level 
staircase N{E) and the averaged integrated density of states N{E) (the smooth dotted line). 

The standard deviation or width cTj of the Gaussians can be taken as a constant for the 
entire spectrum; it is then the free parameter. However, due to the appearance of many 
mini-bands in the Hubbard spectrum for small values of U/t each having different densities 
it is desirable to let CTj adapt to local variations in the spectrum. We have developed the 
following algorithm: take a levels to each side of level i, determine the local average level 
spacing Aj = (-Ej+a — Ei_a)/{2a), and set cr, = 0.608aAj. By this assignment 90% of the 
weight of the broadened peak falls in the interval [E^ — aAj, Ei + aAj] and a becomes the 
free parameter. If a gap (defined as a very atypical spacing) falls within the chosen range 
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we only take the states of the same side of the gap as Ei into account. The procedure is 
illustrated in Fig. |^. We discuss how to optimize the choice of a in Sec. 0. Typically, we 
find a ~ 4 to be a good choice. 

All four methods of unfolding yield essentially the same results. We decided to use the 
Gaussian broadening with varying width, Eq. (||), since it was better suited to the study of 
the Hubbard spectra with its many mini-bands at small coupling strength U/t (see Sec. [V A| ). 

A final remark on the unfolded spectrum is that it is customary to discard from the 
analysis the states closest to the boundary of the spectrum or to the edges of the mini- 
bands. The reason is that these levels in contrast to the levels in the bulk of the spectrum 
do not interact with levels of both higher and lower energy. Hence such levels are nongeneric. 
We usually discarded a few percent of the total number of states on that account. This effect 
is a size effect that is expected to be negligible in the thermodynamic limit. 



B. Quantities characterizing the spectrum 

The simplest quantity one studies in RMT analysis is the probability distribution P{s) 
of unfolded energy spacings s = Ei — Ei-i, where Ei and Ei-i are two consecutive unfolded 
energies. One compares the actual P{s) with the same quantity obtained for random matri- 
ces from one of the matrix ensemble introduced in Sec. |T|. For diagonal random matrices 
P{s) is the Poisson (exponential) distribution P{s) = exp(— s). For NxN GOE matrices 
the distribution of spacings is quite complicated for arbitrary A^, however it is always close 
to the exact spacing distribution of the 2x2 GOE matrices known as the Wigner surmise. 



)GOE/ 



TT 



TT 



— s expl S 

2 ^^4 



(6) 



which therefore is used in practice. 

The spacing distribution probes correlations between consecutive states and is not sensi- 
tive to correlations of higher order. For example an artificial spectrum constructed by adding 
independent variables distributed according to Eq. (|]) will certainly show a Wignerian spac- 
ing distribution but has clearly nothing else in common GOE spectra. To test higher order 
correlations one then looks at the two-point correlation function Y{x) and various weighted 
averages thereof ||^. For GOE in the large A^ limit Y{x) = s{xY + ^^^^ s{t) dt, with 
s{x) = sin(7ra;)/(7rx). One average of Y{x) often studied is the number variance S^(A) de- 
fined as the variance of the number of unfolded energy levels in intervals of length A around 
the unfolded energy Eq: 



NuiEo + ^) - NuiEo - ^) - A 



(7) 



where Nu{e) = I]j^(£ — is the unfolded level staircase, and where the brackets denote 
an averaging over Eq. For the Poissonian case S^(A) = A, while for the GOE case S^(A) = 
A — 2 /o'^(A — x)Y{x) dx with a logarithmic asymptotic behavior. 

Another average of the two-point correlation function is the spectral rigidity A3 (A) de- 
fined as the least square deviation of the unfolded level staircase Nu{e) from the best fitting 
straight line in an interval of length A: 
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A3(A) = l^j min [N^{e) -As- Bfde^ . (8) 

For the Poissonian case the spectral rigidity is A3 (A) = A/15, while for the GOE case 
A3(A) = {\- Jq f{x)Y{x)dx) / 15, with /(x) = {X- xf{2X^ -9Xx - 3x^)/ X'^ and again with 
a logarithmic asymptotic behavior. 

IV. GROUP THEORY AND INVARIANT SUBSPACES 

The problem of diagonalizing the Hubbard Hamiltonian can be reduced considerably 
by group theoretical analysis. Furthermore, as mentioned earlier and as illustrated by an 
example in Sec. |V1 A| it is indispensable for the RMT analysis. The symmetries are explicitly 



dealt with from the beginning of the calculation by constructing the symmetry projection 
operators corresponding to all known symmetries of the model and using them to project 



into symmetry invariant subspaces of the full Hilbert space I^S],^. The main object of this 
Section is to construct the three projection operators Vr, Vs, and Vj corresponding to the 
space symmetry, the spin symmetry and the pseudospin symmetry, respectively. 

A. The space symmetry group 

The first symmetry we consider is the space group Gl of the LxL square lattice with 
periodic boundary conditions. It consists of all permutations g of the sites such that g{i) 
and g{j) are neighbors if and only if i and j are neighbors. In a straightforward manner 
an operator g in the Hilbert space can be associated to each element g of Gl, g\a,b; c, d) = 
\g{a), g{h); g{c), g{d)) , thus forming a group Gl of operators, which commutes with the 
Hubbard Hamiltonian H . For general values of L the space group Gl has been analyzed in 



detail in Ref. |20|. Here we will restrict ourselves to outline this analysis and to correct the 



particular cases of L = 2, which induces a simpler space group, and L = A, which induces a 



much richer space group as briefly mentioned in Ref. . 

For L 7^ 2, 4 the structure of Gl can simply be built up by forming direct |2^ and semi- 
direct 1^ products, denoted ® and (s) respectively, of translation and reflection subgroups. 
Let Tx (Ty) be the subgroups of order L of translations (isomorphic with Zl) in the x (y) 
direction, and let r^, Ty, and be the reflection operations for the x axis, the y axis, and 
the diagonal defined by r^i^x) = —x, ry{y) = —y, and rd{x,y) = {y,x), while e denotes 
the identity transformation. The two subgroups = Tj.(s){e, r^,} and G\ = Tj,(s){e, r^} 
of order 2L (isomorphic with Clv = Zl^Z'i) are formed and combined into the direct 
product subgroup G^l — ® Finally Gl is formed by the semi-direct product Gl = 
G^^(s){e, r^}. One can say that Gl is generated by the elements {tx-,rj.,ty,ry,rfi) although 
this is not the smallest possible set of generators. The order Nl of Gl = {Glv ® Glv)(§)2^2 
is seen to be Nl = 8L^. 

For L = 2 we find G2 = {{e,rri.}^{e,ry})(^{e,rd} = C^^ with = 8. This result differs 
from that of the general case since tx — Tx and ty — Vy. 

For the case L = 4 a richer group appears because the 4x4 lattice with periodic boundary 
conditions is isomorphic with the group of transformations of the four- dimensional unit 
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hypercube. This isomorphism is easily seen by changing the decimal enumeration of Fig. ||a 
into a binary enumeration such that the binary numbers of neighboring sites differs by 
only one bit as shown in Fig. ^ The can immediately be interpreted as the coordinates of 
the corners of the hypercube. Thus the group G4 of neighbor preserving transformations 
are given by combining bit inversion operations 0^1 with permutations of the four bits. 
In Appendix ^ we show that G4 is isomorphic with the group {Z2 ® Z2 ® Z2 ® Z2)@Si 
corresponding to the semi-direct product of all combinations of bit inversions at the four 
positions with the permutation group 5*4 of the four positions. One neighbor-preserving 
transformation of the lattice which is not given by products of the elements {t^, r^, ty, Vy, r^} 
is the hyperplane reflection (see Fig. ^) obtained in 4D by a reflection in the xz-plane 
while keeping y and w fixed. G4 can be generated by replacing above with rh (we note 
that Td = t'^ty^rhtytxVh), and a set of generators is {tx,rx,ty,ry,rh} though in fact as few 
as two elements can be found that generates G4. p8| Instead of 128 elements found by the 
formula = 8L^, we find A^4 = 2^* x4! = 384. We refer the reader to Appendix ^ and 
Appendix p| for further details concerning G4 and Gl, respectively. 




FIG. 3. The mapping of the 4x4 square lattice with periodic boundary conditions onto the unit 
hypercube in 4D. The four double arrows indicate the neighbor preserving transformation r/^. This 
transformation cannot be expressed by the ordinary translations and reflections in 2D, however, it 
can be interpreted as a hyperplane reflection in 4D around the xz plane. 

The complete table of characters Xif \g) of G4 is given in Table |IV| in Appendix while 
those of 6*3, 6*5, and Gg are found by combining Table in Appendix ^ with the method of 
Ref. Armed with the characters tables the character projection operator (where 
subscript R is used to distinguish from S and J) can be constructed for each Zi^/-dimensional 
irreducible representation R': 



(9) 



Using V)i for a given representation together with the Van Vleck basis-function generat- 



ing algorithm P3|,E3] in the actual or a closely related Hilbert space it is straightforward 



to construct an orthonormal //j-dimensional basis {(pi 



(R) 



4"'. 



dR) 



J} and calculate an ex- 
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plicit irreducible representation matrix F-^ [g) = ). Finally the row projection 

operator 'P^\g) can be constructed: 

- 1^ E r[S>-m (10) 

Both projection operators Pjj and Vl^. will be employed in the diagonalization of the 
Hubbard Hamiltonian. 



B. The SU(2) spin symmetry 

It is easily verified that the Hubbard Hamiltonian Eq. (|lD commutes with the z com- 
ponent of the total spin operator Sz = J2fji -S"!*^ , as well as with the corresponding raising 
and lowering operators 5*+ and 5*-. The model therefore possesses a SU(2) spin symmetry. 
We can therefore restrict our diagonalization to the 5*2=0 sector, since the other sectors can 
be reached using the spin raising and lowering operators. The spin operators also commute 
with any space symmetry operator g, so the two symmetry groups form a direct product 
group. The combination of two spin-up and two spin-down electrons under the constraint 
Sz = leads to a total spin quantum number S' = 0, 1 or 2. The explicit construction in 
Appendix y with O = yields the following three spin projection operators Vg^ ^ 

(2|a, b; c, d)+2\c, d; a, b) + \a, c; b, d) + \b, d; a, c) — \a, d; b, c) — \b, c; a, d)) (11) 
{3\a,b]c,d)—3\c,d]a,b)) (12) 
{l\a, b; c, d) + l\ c, d; a, b) — \a, c; b, d) — \b, d; a, c) + \a, d; b, c) + \b, c; a,d)). (13) 

These expressions are generally valid, however, one should note that the Pauli exclusion 
principle reduces the number of terms when electron pairs are present. For example if a = c 
and b = d one finds Vg \a, b; c, d) = Ss'fl\(i, b; c, d). Finally we note, that the spin projectors 
of the state \a, b; c, d) involve all six states generated by permutations of the site indices. In 
Eqs. (|TT|)-(^) only three orthogonal states appear. The remaining three orthogonal states 
are the two S'=l states, |a, c; b, d)—\b, d; a, c) and |a, d; b, c)—\b, c; a, d), and the one S=0 state 
\a, c; b, d) + \b, d; a, c) + \a, d; b, c) + \b, c; a, d). 



vf\aMc,d) = 
Vg'\a, b; c, d) = 
V^s^\a, b] c, d) = 



C. The SU(2) pseudospin symmetry 

The pseudospin symmetry of the Hubbard model has been known for at least a quarter 
of century but recently it was rediscovered and put in the more generalized context of 
the ?7-pairing mechanism [p0| , pl| . As discussed in Sec. || even L lattices can be bi-parted 
using the site sign 9i as an index, and three operators J_, J+, and Jz can be defined: 
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i i \ ia / 

It is seen that J+ creates a pair of electrons with phase 9i on empty sites i. These three op- 
erators form the same algebra as S^, 5*+, and S^, hence the name pseudospin. Furthermore, 
the j-operators commutes with both the space symmetry operators and the spin operators; 



and for the symmetrized Hubbard model H' trivially related to our model if, they 
even commute with the Hamiltonian H' . Hence the (symmetrized) model possesses an ex- 
tra SU(2) symmetry characterized by the quantum numbers J' and Jz analogous to S' and 
Sz for the spin [^]. A detailed analysis shows that the combination of the spin and the 
pseudospin symmetry yields [SU(2)(g)SU(2)]/Z2 = S0(4) rather than the full SU(2)(g)SU(2) 
|3T| . The space symmetry group is completely independent of the spin and pseudospin sym- 



metries. Thus the full symmetry group for L even is ^ = Gl ® SO (4) while for L odd 
it is ^ = Gl® SU(2). In our case the symmetry is lowered in a trivial way from spher- 
ical to cylindrical symmetry in pseudospin space since [?7, J±\ = ±J±, but still we have 
[H, J^] = [H, Jz] = and J' and Jz both remain good quantum numbers. 

For a given lattice with even L containing a fixed number of electrons all states 
have Jz = {Ne — -t/^)/2 (e.g. —6 for L=4 and —16 for L=6). Defining Jq = \Jz\, the 
size J of the pseudospin in this situation takes the values Jo, Jo + 1, • • • , + Ne/2. The 
projection operators Vj^ ^ are found using the construction of Appendix y with O = J^ = 
J+J^ + J^ — Jz- Due to the explicit reference to pairs in the definition of the pseudospin 
operators, it has proven useful to introduce a special notation for one-pair and two-pair basis 
states as follows: 



\P;b,d) = 6'pCp|Cp|4|c)^| I vac) = 9p\P,b; P,d) 



\P;Q) = MQ^^^^^c^rl^ac) = -epeQ\P,Q; P,Q) 
For zero-pair states where a, b, c, and d are all different, only Vj ~ " 



(15) 



Vj^°^\a,b;c,d) = \a,b;c,d), (16) 
for one-pair states where P,b, and d are all different, two projection operators are nonzero, 

V^f°^\P-b,d) = ^^:^^\P-b,d) + -—^ E \Q-^b,d), (17) 
^Wo + ij 2(^Jo + ij Q^p,M 



|P; b,d)= \ \P; b,d)+ \ ^ IQ'Ad), (18) 



while for two-pair states where P ^ Q, all three projection operators are nonzero: 

+ +1V2J +3) ^ ^ 
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2(^Jo + iJl-'o + ^) R^P,QS^P,Q,R 

+ ^7 ^7 ^ E E l^!'^)- (21) 

2 (Jo + 2) (2 Jo + 3) j^^pQg_^pQ j^ 

At this stage all the group theoretical ingredients are ready for the diagonalization of the 
Hubbard Hamiltonian. 



D. The symmetry invariant subspaces 

Using the projection operators the A'^iif-dimensional Hilbert space can be broken down 
into smaller symmetry invariant subspaces. Let Q be the group containing the symmetry 
operations 7, each being a product of a space symmetry transformation, a spin rotation and 
(for L even) a pseudospin rotation, 7 = g®gs®9j- Q thus consists of one finite group and one 
or two compact Lie groups. Let p be a multi-index describing an irreducible representation 
r('') of Q. For even (odd) L we have p = {R',S',J') {p = (R'^S')). The "celebrated" 
theorem [^,0 states how many times Op each row of the irreducible representation T^^^ 
with character x^^^ appears in any given not necessarily irreducible representation F with 
character x- 



I d^X^P>\^)x{l) = ^J2[ dgs f dgjx^p>{^)x{l) (22) 
JG Nr. rr^ Jsv(2) Jsv(2) 



Jg 1\l JSV(2) JSU(2) 

We now chose F to be the following A'^n-dimensional reducible representation: 



F,,(7) = (^|7|J) x{l)=T.^^m■ (23) 



1=1 

Inserting into Eq. ( p2|) these x(7)'s which are directly linked to the Hilbert space yields: 

Nh a Nh 1 Nh 

„ - . f E(^I^^^^K) = f „ 



Nh 1 Nh i Nh 

I ^7x^^^*(7)E(^l7K) = r E(^I^^'^K) = r E(^l^f ^ ® ^f^K), (24) 

JQ .-I tR ._i fR .-I 



which is obtained by using ^(7) = x{.g®gs®gj) = x{g)x{gs)x{gj) the definitions 
of the character projection operators. Note that only the space group projection needs a 
normalizing factor 1/Ir. The spin and pseudospin are both restricted to take only one value 
of their respective z component; hence their normalizing factor is 1. 

The initial size of the Hubbard Hamiltonian to be diagonalized is [L(L — 1)/2]^. However, 
only states transforming according to the same row of the same irreducible representation 
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p can have a nonzero matrix element. Hence the Hamiltonian matrix breaks up in blocks, 
one per representation, and within each representation a further division into Ir equivalent 
blocks occurs. The relationship between Nh, cl{r',s',j'), and is = J2{r,s,j) ^r'0'{R',s',j')- 
This reduction is considerable as shown in Table |. 





L = 3 


L = 4 


L = 5 


L = Q 


Nh 
maxp{ap} 
Ep^l/Nfj 


1296 
38 

14.1x10"^ 


14400 
146 
0.7x10"^ 


90000 
1794 
2.3x10"^ 


396900 
5490 
1.0x10"^ 



TABLE I. For L = 3, 4, 5, and 6 are shown the dimension Nh of the total unreduced Hilbert 
space and the dimension maxjop} of the largest symmetry invariant subspace. Since matrix diag- 
onalization is a n^-operation, the number '^pa^/Nfj provides an estimate of the relative reduction 
in computer time by projecting into the invariant subspaces. 



V. NUMERICAL DIAGONALIZATION OF THE HUBBARD MODEL 

The numerical calculation of the exact spectra of the Hubbard model begins by determin- 
ing the block size Op for all the irreducible representations p. Then for a given p = {R', S', J'), 
where J' is disregarded in case of odd L, sua ap- dimensional orthonormal basis is found by 
using the projection operator 

= ® vP ® \ (25) 

which projects onto the first row of the space group representation R' with the spin index 
S' and pseudospin index J'. The projection operator Vq is applied on one basis state after 
another while performing an ongoing Gram-Schmidt orthonormalization procedure. This 
yields new basis states \wi), and when ap such states are found, the procedure is terminated. 

Next step is to calculate the kinetic energy and potential energy matrix elements 
{wi\T\wj) and {wi\U\wj). These matrix elements are stored in the computer and thereafter 
it is a simple matter to pick any value of U/t and diagonalize the Hubbard Hamiltonian, 
H = —tT + UU , using standard diagonalization routines. The calculation of the symmetry 
invariant matrix blocks of T and U takes time of the order of one diagonalization. 

A. The spectrum and first order perturbation theory of the ground state 

In this section we discuss some features of the raw spectrum of the Hubbard model. Only 
later in Sec. [V^ we are going to unfold the spectrum and look for universal features. As a 
function of the coupling strength U / 1 the spectrum clearly falls in three classes. In the weak 
coupling limit {U/t <^ W/t), to be studied in more detail below, the spectrum acquires a 
band width W ~ 32t. It consists of a number M of well separated mini bands reminiscent 
of the huge degeneracy at zero coupling due to size quantization. For intermediate coupling 
strengths {U/t ^ W/t) the spectrum becomes rather featureless. No apparent gaps or bands 
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emerge. In the strong coupling limit {U/t ^ W/t) the spectrum splits up in three well 
separated bands centered around the energies E/t = 0,U/t,2U/t. These are the Hubbard 
bands corresponding to states containing zero, one or two pairs of electrons. In Fig. |^ are 
shown one typical spectrum from each of the three coupling strength regimes. 




-10 -5 5 10 -10 -5 5 10 15 20 500 1000 1500 2000 

E/t E/t E/t 



FIG. 4. The integrated density of states N{E) of the Hubbard model with 4 electrons on 
the 5x5 lattice for the irreducible representation {R, S) = (10,0) with U/t = 0.1,10,1000. For 
U/t = 0.1 the mini bands are still clearly visible, N{E) is essentially featureless for U/t = 10, 
while for U/t = 1000 the three Hubbard bands at 0, 1000 and 2000 appears. The inserts are 
magnifications of N{E) where the arrows point. The smooth N{E) is added. 

We now focus on the weak coupling limit where it is natural to make Fourier transforms 
from real space to momentum space, c|^^ = X]j exp(zk- rj)c]^, where k = {k^,ky), with 
/c^ = 27m /L and n = 0, 1,...,L — 1. The eigenstates of the kinetic energy operator —tT are 

|ko,ki;k2,k3) = cl^^cl^^cl^^cljvac) 

3 

-tf |ko, ki; k2, kg) = -2t J2 [cosiK) + cos(A:^)] |ko, ki; ks, kg). (26) 

n=0 

The Pauli principle prevents ko=ki and k2=k3, and the ground state energy e][^^ becomes 

= -12t - At cos{2tt/L). (27) 

For large lattices this tends toward —16t and results in a band width W = 32t. For L = 
3, 4, 5, and 6, respectively, the actual band widths W/t are 20.0, 24.0, 26.2, and 28.0 while 
the number M of mini bands for U = are 7, 13, 42, and 29. 

The limitation of perturbation theory is demonstrated by a first order degenerate per- 
turbation calculation for the ground state. Let pi be the one-dimensional momentum com- 
ponent pl = 2tt/L and construct the five two-dimensional momentum vectors = (0,0), 
q° = (pl,0), = {0,Pl), = (— Pl,0), and q^ = (0,— Pl). The 16 states \n,X), defined as 

l/i. A) = |0t, q^; 0^, qj), /i, A = 0, 1, 2, 3, (28) 

form the ground state multiplet separated from the next multiplet by an energy gap AEl = 
2t[l — cos(pl)]. In momentum space the interaction operator ULf takes the form 

UU=^ 4i+qTCkiTCk2-qiCk2i, (29) 

ki,k2,q 
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and after some simple algebra the matrix elements of the perturbation are found to be 

(m, l\UU\fi, ^) = J2 (3'^m,M'^Z,A + ^m+l,^L+\) ■ (30) 

The eigenvalues can be found analytically, and we end up with the following expression for 
the perturbed levels Ei^^piJJ) with degeneracies d: 

El^0{U) = Ef* + (3—, /3 = 3(rf=7),4(rf=4),5(rf=4),7(rf=i). (31) 

Thus in first order perturbation theory the ground state degeneracy is partly lifted leaving 
a 7-fold degenerated ground state for [/ > 0. 

In Fig. ^ we compare the exact numerical calculation with Eq. ( |3T| ) for L = 5 and 6. 
Note especially how in the exact calculation the degeneracy of the ground state energy is 
lifted completely. This is also demonstrated in Table |T|, where it can be seen that also for 
L = 3 the ground state is nondegenerate. However, the ground state of L = 4 remains 
3-fold degenerated even in the exact calculation. This reflects the particular symmetry of 
the 4x4 lattice, where the hyperplane reflection leads to the existence of 3-fold degenerate 
representations as shown in Appendix |^. 




0.00 0.02 0.04 ,0.06 0.08 0.10 



FIG. 5. The evolution of the 16-fold degenerated ground state multiplet as a function of U jl? . 
The exact numerical results for L = 5 (dotted lines) and L = 6 (dashed lines) are contrasted with 
the first order perturbation calculation (full straight lines). The insert is a magnification of that 
portion of the plot marked by a rectangle showing how the exact calculation leads to a splitting into 
three sub-multiplets with degeneracies 1, 2, and 4 of the 7- fold degenerated perturbation theory 
ground state. The exact ground state is nondegenerate. See also Table 

Using the quantum numbers {k^^ ky, b^, by, c) of Table |V| in Appendix ^ to interpret the 
representation label R in Table || we note that they are of the same form for L = 3,5, 6: 
(0,0,1,1,-1), (0,0,1,-1,*), (f ,f ,*,*,-l), (f ,0,*,1,*), (f ,f ,*,*,1), and (0,0,1,1,1), where the 
first set corresponds to the unique ground state. This result extends to the L=4 lattice 
when the actual space group is replaced by the one generated by {tx,rx,ty,ry,rd}. One 
state in the L=4 triplet ground state belongs to the (0,0,1,1,-1) representation and the other 
two states to (^,0,1,1,*). In all cases the ground state is odd with respect to the diagonal 
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reflection r^, which is understandable since that symmetry suppresses the abihty of having 
pairs along the diagonal. Finally we note that due to the aforementioned spin symmetry of 
the Hubbard model, the ground state of the 5*^=0 sector is in fact a global ground state, 
and furthermore since it has S = the same energy level does not exist in any other 
sector. Thus we can conclude that the global ground state in the generic case (L 7^ 4) is 
nondegenerate. 





L = 3 


L = 4 


L = 5 


L = 6 


# 


E/t {R,S) Ir 


E/t ( R,S,J) Ir 


E/t ( R,S) Ir 


E/t ( R,S, J) Ir 


p 


-9.94000 ( -, -) 7 


-11.94000 ( -, -, -) 7 


-13.17607 ( -, -) 7 


-13.94000 ( -, -, -) 7 




1 

2 
3 
4 
5 


-9.94122 ( 2,0) 1 
-9.94110 (8,1) 4 
-9.94097 (4,1) 2 
-9.92188 ( 5,0) 4 
-9.90229 ( 7, 0) 4 
-9.86204 ( 0,0) 1 


-11.94268 ( 9,0,6) 3 
-11.94232 (16,1,6) 6 
-11.90425 (17,0,6) 6 
-11.86560 ( 0,0,6) 1 


-13.18059 ( 2,0) 1 
-13.18012 ( 4, 1) 2 
-13.18005 (10, 1) 4 
-13.16204 ( 7,0) 4 
-13.14302 ( 9,0) 4 
-13.10590 ( 0,0) 1 


-13.94676 ( 4,0,16) 1 
-13.94617 ( 9,1,16) 2 
-13.94606 (25,1,16) 4 
-13.92870 (14,0,16) 4 
-13.91032 (24,0,16) 4 
-13.87494 ( 0,0,16) 1 



TABLE II. For U/tL"^ = 0.02 the 7- fold degenerated perturbation theory ground state i?}^ 
given in row P is compared with the exact energies E of the 16 levels splitting off from the C/ = 
ground state E^j^\ The degeneracies of the sub-multiplets are given by the dimension Ir of the 
corresponding irreducible representation R of the space group. Note that except for L=4 the exact 
ground state is nondegenerate. Also listed are the quantum numbers (i?, S, J) of the levels. 



B. Remaining degeneracies 

After taking all the known symmetries into account and after projecting into the sym- 
metry invariant subspaces, it turns out that within each subspace some further degeneracies 
remain. The reason is that the low fllling of the lattice allows for a kind of restricted per- 
mutation symmetry for the particle momentum components resulting in energy eigenstates 
which are simultaneously eigenstates of T and U and therefore independent of U (though 
their eigenvalues might depend on U). We denote such states T /U states or iV's), where S is 
the spin and 7 the eigenvalue of tJ . With four electrons 7 takes the values 7 = 0, 1, 2 corre- 
sponding to superposition of states containing exactly 7 pairs, i.e. doubly occupied sites. A 
generic energy eigenstate is not a T/f/ state and in that case we write 7 = *. In Table |ITT| we 
show the number of generic energy eigenstates and T /JJ states found numerically. The Hub- 
bard model is thus partly integrable. The spectrum of each symmetry invariant subspace 
is a mixture of an integrable component (the T /U states) and a non- integrable component 
(the generic states), and in the analysis in Sec. of the spectral statistics we throw away 
the integrable component and analyze only the generic non-integrable component. 

The T/f/ states are formed by specific superpositions of eigenstates of T by permuting 
the eight momentum components A;^, ( = x,y and n = 0, 1,2, and 3 such that the sum of 
cosines in Eq. (^) remains unchanged. In the following we mention some large classes of 
such states. The reader is referred to Appendix for details. 
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7 


S 


L = 3 


L = 4 


L = 5 


L = 6 







540/(0) 


4169/(0) 


31977/(0) 


115896/(0) 


* 


1 


621/(0) 


4472/(0) 


41662/(0) 


137199/ (0) 




2 


0/(0) 


0/(0) 


0/(0) 


0/(0) 







0/(0) 


1176/(1143) 


523/(316) 


23555/(23427) 





1 


n / /r\\ 

9/(0) 


2548/(2519) 


3188/(2823) 


60306/(60160) 




2 


126/(80) 


1820/(1727) 


12650/(12362) 


58905/(58723) 







0/(0) 


94/(14) 


0/(0) 


408/(247) 


1 


1 


0/(0) 


120/(24) 


0/(0) 


630/(420) 




2 


0/(0) 


0/(0) 


0/(0) 


0/(0) 







0/(0) 


1/(0) 


0/(0) 


1/(0) 


2 


1 


0/(0) 


0/(0) 


0/(0) 


0/(0) 




2 


0/(0) 


0/(0) 


0/(0) 


0/(0) 




% 


89.6% 


60.0% 


81.8% 


63.8% 



TABLE III. For L = 3,4, 5, and 6 are shown the number of energy eigenstates found nu- 
merically in each of the four main groups indexed by 7. The first group, denoted 7=*, contains 
the generic energy eigenstates that are not eigenstates of tJ . The three other groups, denoted 
7 = 0, 1,2 respectively, contain the nongeneric T /tJ states that are simultaneous eigenstates of U 
(with eigenvalue 7) and T. In parenthesis are given the number of states that remain degenerated 
in the symmetry invariant subspaces after taking the space, spin, and pseudospin symmetries into 
account. Note how the 7=* states exhibit no further degeneracy. In the last row denoted "%" are 
given the percentages of 7=* states out of the total number of states. 

First we note that for even L proportionally many more states remain degenerate as 
compared to L odd. This difference is due to the momentum component k'^ = vr, which only 
exists for L even. Because cos(7r — A;) + cos(A;) = independent of k such terms drop out 
when T is applied to a state containing this combination and a certain degree of freedom 
is left to form the T/U state superpositions. Starting with two-pair states 7 = 2 we find 
exactly one energy eigenstate {ip^Z^). It depends on the momentum vr, and hence it exits 
only for even L. Similarly, for one-pair states only even L leads to T/U states. A number 

{\ ) of states \ipsZi) can easily be constructed. More care must be taken upon forming 
the corresponding S=0 states. However, we have succeeded in constructing analytically the 
number of states \ipsZo) required by Table |IT1[ Finally, for zero-pair states the existence 
of the momentum tt is not required to form the T/U states, but it certainly helps. Thus 
both even and odd values of L lead to nongeneric energy eigenstates, but relatively more 
such states are found for even L. It is easy to see that all states with the maximal spin 
S = 2 are T/U states {ip'^'I^). This is a trivial consequence of choosing four different sites 
(or momenta) and forming the superposition given in Eq. (|T3]). There are (^ ) such states. 
The construction of states lifjgZ^g,) with S'=l or is more cumbersome, so in Appendix 
we have only given two examples of classes of such states. 
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The main result of this section is that the degeneracies that remain after space, spin and 
pseudospin symmetry have been taking into account, is related to a restricted permutation 
symmetry of the momentum components. We have not found the associated projection 
operators, but numerically and partly analytically we have established the fact that the 
degenerate states are simultaneously eigenstates of T and U with energies E{U) = E{0)+jU. 
By discarding these states we end up with nondegenerate ^/-dependent states. It is the 
spectral statistics of these states we analyze in the following section, and this analysis confirm 
our claim that all symmetries indeed have been taken into account. 



VI. SPECTRAL STATISTICS OF THE HUBBARD MODEL 



Having sorted the spectrum according to all symmetries including the restricted permu- 
tation symmetry of the momentum components the RMT analysis can be performed. As 
discussed in Sec. |ITT| the first step is unfolding of the spectrum. There we mentioned how 
the unfolding procedure is not uniquely determined, so we turn to that problem first. 



A. Optimization of the unfolding procedure 




0.0 

s X 
FIG. 6. The spectral statistics P{s), A3(A), and S^(A) for the invariant subspace {R, S) = 

(13,1) of the 5x5 lattice as a function of the free parameter a of the unfolding procedure using 

Gaussian broadening with variable width. In all three panels the values for a are 2, 4, 5, 6, 8, 10 

and 20. In the first panel 6 of the 7 curves fluctuates around P^^^[s) hardly visible as a smooth 

dotted line. In the two other panels the Poisson case are given as a dotted straight line while the 

GOE case is given by a dotted smooth curve. For all data sets a few representative error bars are 

shown. 



We unfold the spectra by using the method of Gaussian broadening with variable width 
discussed in Sec. [Ill A| . The free parameter a corresponding to how many energy levels each 
Gaussian essentially spreads out over to each side, and the question now arises how to choose 
it. The problem we face is illustrated in Fig. ^ were it is seen that although the level spacing 
distribution P{s) is essentially independent of the choice of a, both the number variance 
E^(A) and the spectral rigidity A3(A) varies with a. It is seen that A3(A) is less sensitive 
to changes in a than S^(A). The former seems to saturate for large values of a, while the 
latter continues to grow rapidly as a enhances. We have chosen that value of which 
makes S^(A) fit the corresponding GOE curve as well as it can. In general that leads to 
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ao ~ 4. We note that for this choice of a, A3 (A) is almost saturated, while P{s) remains 
unchanged. Thus, two of the three statistics are essentially independent of variations of a 
around a^, while the third is as close to the GOE behavior as it can be. 

To illustrate the importance of sorting the spectrum by group theory we show in Fig. |^ 
the level spacing distribution for L = 4 with U/t = 10 for the case where all symmetries are 
taken into account (full symm.) and for the case with lower symmetry (low symm.) where 
the spin and pseudospin symmetries are being kept intact but where the space group has 
been reduced by replacing among the generators the special hyperplane reflection with the 
ordinary diagonal reflection r^. For the full symmetry case a distribution rather close to the 
Wigner surmise is found, whereas the level repulsion is partially lost for the low symmetry 
case, and the data fits reasonably well the distribution found by mixing two GOE-spectra 



2^ with relative weights 0.72 and 0.28. This makes sense since by lowering the symmetry 



group artificially, spectra from the independent true symmetry invariant subspaces are being 
mixed. For the more severe symmetry reduction where the pseudospin is altogether neglected 
and the space group is generated only by {tx,ty}, we find the level spacing distribution to 
be the Poisson (exponential) distribution (not shown in the figure). 




s 

FIG. 7. P{s) for L = 4 with U/t = 10 is calculated after sorting the spectra using either the full 
symmetry group of the Hamiltonian (o) or an symmetry group artificially lowered (*) as described 
in the text. The full symmetry case compares well with the Wigner surmise (smooth full curve) 
whereas the lower symmetry case compares well with the distribution of two GOE-spectra mixed 
with the relative weights 0.72 and 0.28 (dotted smooth curve). 



B. The statistics P{S), ^^(A), and A3(A) of the Hubbard model 



In this subsection we present the results of the spectral statistical analysis of the Hubbard 
model at low filling. We present only results for L=4, 5, and 6 since L=3 yields too poor 
statistics due to its small invariant subspaces. Besides letting L vary we are also varying 
the coupling strength U/t and present results for weak, intermediate and strong coupling 
regimes for L = 5. To improve on the statistics we have averaged over the largest invariant 
subspaces. For P{s) the size SP of the error bars shown in the figures are estimated by 
6Pi = C^/rii/hi, where rii is the number of points in bin i of the associated histogram, hi 
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is the width of the bin, while C is the normahzation factor rendering a total probability of 
1. For S^(A) and A3(A) the error bars are estimated by ordinary variances obtained from 
the values calculated in the many contiguous intervals of length A throughout the unfolded 
spectrum. 

The results for P{s) are shown in Fig. ^. It is seen that for all lattice sizes and for any 
value of the coupling strength the level spacing distribution is fairly close to the Wigner 
distribution of GOE; it possesses a pronounced linear level repulsion for small s, a peak near 
s = 0.8 signaling spectral rigidity, and a rapid fall off for s > 2. 




1 2 3 1 2 3 



s s 
FIG. 8. The probability distribution P{s) of the level spacings s averaged over the largest 

symmetry invariant subspaces for the Hubbard model with four electrons on square LxL lattices. 

To the left is shown P{s) as a function of lattice size L at medium coupling strength. The data 

represents L=4, 5, and 6 for U/t = 10.0. To the right is shown P{s) as a function of coupling 

strength at fixed low filling. The data represents U/t=0.1, 10, and 1000 for L=5. The full line is 

the Wigner distribution found for GOE random matrices. 




02468 10 02468 10 

X X 

FIG. 9. The number variance S^(A) calculated for the same parameters as in Fig. |^. The data 
are compared to the results of the random diagonal matrix ensemble (Poisson) and the random 
full matrix ensemble (GOE), shown as the straight full line and the curved full line, respectively. 

In Fig. ^ is shown S^(A) of the Hubbard model for the same parameters as for P{s) just 
mentioned. When A is small the rigidity of the Hubbard spectrum is very close to that of 
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the GOE random matrices, while for larger A a saturation sets in. For all values of U/t we 
find the critical value A* where the departure from GOE sets in to be roughly 2. The precise 
origin of A* remains unclear. 

Finally, in Fig. |10| are shown the results for the spectral rigidity A3 (A). As S^(A) also 
A3(A) displays an excellent agreement with GOE for A < A* ~ 2, for all fillings and for 
all values of U/t. The deviations from GOE beyond A* are not so marked. The curves lie 
between the Poisson line and the GOE curve, but rather close to the latter. 
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02468 10 02468 10 

FIG. 10. The spectral rigidity A3 (A) calculated for the same parameters as in Fig. ^. The data 
are compared to the results of the random diagonal matrix ensemble (Poisson) and the random 
full matrix ensemble (GOE), shown as the straight full line and the curved full line, respectively. 



It is remarkable how the results for the three statistics studied are fairly independent 
of the size of the lattice (equivalent to the filling) and of the coupling strength. We find 
GOE-like behavior not only for all finite values of U / 1 including those close to the integrable 
f/ = limit, but also for filling factors as low as 0.06 close to the integrable single-particle 
limit. However, as is evident from the behavior of especially S^(A) at large energy scales, 
the Hubbard model cannot be modeled exactly by a simple GOE random matrix model. 



VII. CONCLUSIONS AND DISCUSSION 

In this paper the energy level statistics of the Hubbard model for L x L square lattices 
[L = 3,4,5,6) at low filling (four electrons) has been studied numericaUy for a wide range 
of the coupling strength. With great care all known symmetries of the model (space, spin 
and pseudospin symmetry) have been taken into account explicitly from the beginning of 
the calculation by projecting into symmetry invariant subspaces. The details of this group 
theoretical treatment was presented with special attention to the nongeneric case of L = 
4, where a particular complicated space group appears. The resulting reduction of the 
numerical diagonalization is significant, and the method presented can in a straightforward 
manner be extended to larger lattices and higher fillings and thus form the basis of improved 
numerical studies of the Hubbard model and related models without disorder. In particular, 
this work can be used as a starting point for calculating various spectral functions, for which 
explicit forms of the eigenstates are required. This will be dealt with in forthcoming work. 
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For all the lattices studied a significant amount of levels within each symmetry invariant 
subspace remain degenerated, but except for L = 4 the ground state is nondegenerate. We 
explained the degenerate states as a consequence of a restricted permutation symmetry of 
the momentum components. These states, all independent of U, form an integrable part of 
the spectrum, and after discarding them we end up with nondegenerate spectra on which 
the level statistical analysis could be performed. 

The intricate structure of the Hubbard spectra necessitated the development of a careful 
unfolding procedure as a preparatory step before the level statistical analysis. The procedure 
we arrived at tested favorably in many cases of pathological spectra, and it seems to be very 
robust and applicable in general cases were no other natural unfolding procedure exist. 

Finally, we have performed a level statistical analysis of the Hubbard spectra, and we 
presented results for the level spacing distribution P{s), the number variance S^(A), and 
the spectral rigidity A3 (A). The statistics for the different lattice sizes and for a wide range 
of coupling strengths are essentially the same: P{s) shows a good agreement with GOE. 
S^(A) agrees only with GOE up to the ?7/t-independent medium sized energy scale A* ~ 2 
beyond which a saturation sets in. A3 (A) also agrees with GOE for A < A*. The deviation 
from GOE beyond A* is not so marked as that for S^(A). The curve falls between that of 
GOE and Poisson but rather close to the former. We stress that these results were also 
obtained for very small coupling strengths approaching the integrable zero coupling limit. 
This emphasizes the nonperturbative nature of the model revealed by our analysis: even the 
smallest deviation from the integrable limits leads to spectral statistics usually associated 
with non-integrability and quantum chaos, and in this sense the model seems always to be 
in the strong coupling limit. 

Largely, our results show GOE-behavior of the spectral statistics of the typical high lying 
excitations of the Hubbard model at low filling. This indicates that at least the incoherent 
part of the electronic spectral functions (related to the coherent part describing the low 
lying electronic excitations through sum rules) are out of reach by standard methods. On 
the other hand, it should be possible to model this part by a random matrix ansatz. This is in 



agreement with previous results of the tJ model near half filling |T^. However, the cause of 
the deviations from GOE we found in S^(A) and A3 (A) beyond A* remains an open question. 
The similar question has been answered in general for single-particle systems with mean-level 
spacing A: in disordered (metallic) systems A* A ~ h/ro-, where To is the time it takes a 
particle to diffuse through the system whereas for pure (ballistic) systems A* A ~ ^/tq, 
where Tq is the period of the shortest periodic orbit . Results are also beginning to emerge 
for disordered interacting lattice systems, where A* is related to the ratio U /W between some 
interaction strength U and the disorder induced single-particle band width W [TB,3B|. For 
these systems the deviations from GOE are due to the preferential basis supplied by the 
given disorder potential. For example the system consisting of two interacting particles 
in a disorder potential [^] could be studied analytically by adding a random diagonal 
matrix modeling the disorder single-particle states to a random GOE matrix modeling the 
interactions between these states, and it was found that S^(A) increased as a power law for 
A > A*. This behavior is in contrast to the saturation we found (see Fig. which looks 



more like the result of the ballistic single-particle case [^]. It is perhaps not surprising 
that such a similarity exists between the disorder-free chaotic single-particle case and the 
disorder-free Hubbard model rather than between two strongly correlated systems one with 
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and the other without disorder. However, exactly what physical mechanism produces a 
preferential basis for the Hubbard model at low filling is not known, and neither is it known 
why a similar mechanism is suppressed for the tJ model near half filling, where much less 
pronounced deviations from GOE are found |]15|- These questions are topics for future work. 



VIII. ACKNOWLEDGEMENTS 

It is a pleasure to thank Benoit Dougot, Jean-Louis Pichard, Clement Sire and Dietmar 
Weinmann for stimulating discussions and Elliot H. Lieb for valuable comments. We thank 
the Centre Grenoblois de Calcul Vectoriel for its kind help with the part of our numerical 
work that was carried out at its CRAY-94 computer. H.B. was supported by the European 
Commission under grant no. ERBFMBICT 950414. 



APPENDIX A: THE IRREDUCIBLE REPRESENTATIONS OF Ga 



In Sec. |1V A| we showed how the group G4 of neighbor-conserving transformation of the 



4x4 lattice is isomorphic with the point group of the four-dimensional hypercube. In this 
appendix we determine the structure of this group and we sketch how the irreducible rep- 
resentations are found analytically. The fundamental simplification is the observation that 
6*4 has the structure of a semi-direct product involving an invariant Abelian subgroup. The 
theorem of induced representations ||2^ can then be used to find all irreducible representa- 



tion. The theorem is stated below. The coordinate set of a corner in the unit hypercube in 
four dimensions is given as a four bit binary number. Any transformation of the hypercube 
can thus be written as a permutation of the four bits followed by bit-inversion (0 1) of 
all, some or none of the bits. 

In what follows any quadruple (xi, X2, ^3, 0:4) is written in short hand notation as (xj). 
The group of bit permutations is of course the permutation group 5*4 denoted B in the 



following to be consistent with Ref . , on which the group theoretical work in this appendix 
is based. Any element h&B is written as 6 = (pi) = (61, 62, ^3, ^4)5 listing the permutation 
of the numbers 1, 2, 3, and 4. The group of bit inversions is denoted A. It is easily seen 
that ^ = ^2 ® ^2 ® ^2 ® 2^2, since bit- inversions does or does not take place on each of the 
four bit positions. Any element a&A has the form a = (a^) = (ai, a2, a^, 04), with Oj = ±1, 
where means no bit-inversion and —1 means bit-inversion. Any element of g can 
be written as g = ah = (aibi) and conversely all products ab^G^. A contains 16 elements 
and B contains 24 so that contains 384 elements. It is readily verified that A is an 
Abelian subgroup of 6*4. Furthermore, A is an invariant subgroup since for WaEA,WbEB: 
bab~^ = b{ai)b~^ = {ab(i))bb~^ = (afe{i)) eA. Finally, the only common element of A and B 
is the identity. We can therefore conclude that 6*4 is a semi-direct product of the invariant 
Abelian subgroup A with B: 

G4 = A(§)B = {Z2(g)Z2(S)Z2(S)Z2)(§)S4. (Al) 

The first step in calculating the irreducible representations for G4 is to construct the 
character table of A. This is easily found as the product of the character table n~J" 
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for Z2 with itself four times. The 16 irreducible representations are identified by the index 
q = (f?!, 92, (hi Qi), with gj = or 1, and the q-th character x'^{(^) is given by the product of 
to the power qf 



1=1 



Next step is to pick any q' and construct the associated little group i3(q') C B defined as 

i3(q') = {beB\ x'^ipah-^) =x'^{ci), for VaG^}. (A3) 
The group B is then written in a coset decomposition after i3(q'): 

B = B{c{)hi ® i3(q')62 © ... © S(q')&M'. (A4) 
For each of the M' coset representatives hj an index q'^ is determined such that 

X"^^ = x'^'ibjabj^), for VaG A (A5) 

Note that bi is the identity and that q\ hence equals q'. The set orb(q') = {q'l, q'2, . . . , qV'} 
is called the orbit of q'. Now pick a q" outside orb(q') and repeat the procedure. This is 
continued until each of the 16 q's are associated with an orbit. The last thing to do before 
constructing the irreducible representations of G4 is to find the irreducible representations 
A'l P of the little group i3(q'). This is usually a simple step due to the small size of B{q'). 
The theorem of induced representations states that all irreducible representations 



V^P of the semi-direct product A@B, A being an invariant Abelian subgroup, are found 
as follows, (i) Pick one q' from each orbit, (ii) Construct of the little group i3(q') and its 
Tip irreducible representations A'^P,p = l,...,np. (iii) Find the M' coset representatives 
bj GB,j = 1, . . . , M' of B with respect to -B(q'). (iv) Then the matrix elements of T^'^ for 
element ab is given by: 

^ ' lo iibkbbf ^B{cl). ^ ' 

Here we will not give the explicit expressions of the irreducible representations of 6*4, 
but rather just briefly sketch the construction of them and calculate how many there are 
and what is the dimension of each of them. 

First we chose q' = (0000). From Eq. ( [A^ ) it is easily seen that (a) = 1 for 

\/aeA. Hence x^^^^^\bab-^) = for Va G A,ybeB, and according to Eq. (P) we 

find S(OOOO) = B. The coset representation of B consists of only one term and the single 
coset representative 61 is the identity. As a consequence we have orb(OOOO) = {(0000)}. 
Finally, the irreducible representations of i3(0000) are simply those of ;B = 6*4 (or 
as the group also is called |2^), i.e. there are 5 different A^^-matrices with dimensions 



1, 1, 2, 3, and 3 respectively. From Eq. ( [A 61 ) we find that j,k = 1 and thus we have 
found 5 irreducible representations of G4 with dimensions 1, 1, 2, 3, and 3 of the form 
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TABLE IV. The character table of G4 defining the index R for each of the 20 irreducible 
representations and the index C for each of the 20 classes. Column C = contains the dimension 
Ir ranging from 1 to 8 of the representations. 

Next we choose q' = (1111). We find that for Va G A,ybeB: x^^^^^\bab-^) = = 
Ylj ttj = x^^^^^Ko.)- So in analogy with q' = (0000) we have ;B(1111) = B, and like before 
we find 5 irreducible representations of G4 with dimensions 1, 1, 2, 3, and 3 of the form 
_ -1^(1111) ((2)/\0p^^^^ with the same A°^'-matrices but different x-prefactors as for 
q' = (0000). 

We go on with q' = (1000). This yields x^^^^^Kbab'^) = a6(i) which equals x^^°°°H«) 
for Va G ^ if and only if 6(1) = 1. Thus i3(1000) is the six-element subgroup of B which 
leaves the first axis invariant. The coset decomposition of B contains four terms with coset 
representatives that each leaves one of the four axis invariant. Not surprisingly we find 
orb(lOOO) = {(1000), (0100), (0010), (0001)}. The little group ;B(1000) is isomorphic with 
C31, and has thus 3 irreducible representations A^'^' with dimensions 1, 1, and 2 respectively. 
This combined with the fact that index j, k in Eq. { \KB^ ) runs over the four coset repre- 
sentatives means that we have found 3 more irreducible representations r*^^°°'')^ of G4 with 
dimensions 4, 4, and 8, and with entries of the form 0, ;;^(iooo)^ip^ ^(oioo)^ip^ ^(ooio)^ip^ 

^(OOOl)^lp^ 

Then we consider q' = (0111). The character x^' now gives x^^^^^Hbab~^) = ab(2)ab{3)ab{i) 
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which as above equals x^^^^^\(^) for VaG^ if and only if 6(1) = 1. The rest of the analysis 
is similar as the previous case: orb(Olll) = {(0111), (1011), (1101), (1110)} and S(IOOO) is 
isomorphic with C^^. We end up with 3 more irreducible representations r'^^^^^-'^ of G4 with 
dimensions 4, 4, and 8, having entries of 0, ;^(oiii)a1p, ;^(ioii)/^ip^ ^{iioi)^ip^ ^{iiio)^ip_ 
The last choice for q' turns out to be q' = (0011), and we obtain x^'^^^^Kbab~^) = ci6(3)afe(4), 
which equals x'-°°^^-'(a) for Va G ^ if and only if b does not mix the pairs (1,2) with (3,4). 
It is easily seen that ;B(0011) is a 4-element group isomorphic with Z2 ® and thus has 
4 1-dimensional irreducible representations A^^. The coset representation of B consists of 6 
terms in this case, and orb(OOll) = {(0011), (0101), (0110), (1001), (1010), (1100)}. At this 
point we note that all 16 possible values of q' now is a member of an orbit. The 6 coset 
representatives combined with the 4 1-dimensional irreducible representations of ;B(0011) 
yields through Eq. ( |X6| ) 4 new irreducible representations of 6*4 each being 6-dimensional and 
each having entries of the form 0, ;^(iioo)^2p^ ^(ioio)^2p^ ^(iooi)^2p^ ^(oiio)^2p^ ^(oioi)^2p^ 

or t^(o°i^)A^''. In conclusion we see that G4 has 20 irreducible representations with the 
dimensions listed in the first column of the character table shown in Table ITVI. This result 



is to be contrasted with the list of representation dimensions given in Ref. ||20|, where only 
translations and reflections are taken into account in the analysis of the space group of the 
4x4 lattice. Note especially the 3- and 6-dimensional representations found here as opposed 



to dimensions equal powers of 2 in Ref. |20]. 



APPENDIX B: THE IRREDUCIBLE REPRESENTATIONS OF G3, G5, AND Ge 
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The 5x5 lattice 
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The 6x6 lattice 
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TABLE V. The irreducible representations R of Gl for L = 3, 5, and 6, their corresponding 
quantum numbers {kx,ky,bx,by,c), and their dimensions Ir. The symbol * refers to indefinite 
reflection quantum numbers. 
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The irreducible representations R of the space groups Gl with L = 3, 5, 6 can be derived 
analytically following Ref. |20|. In Table ^ each of them are specified by the lattice size 
L and by their representation quantum numbers {kx,ky,bx,by,c) related to the translation 
and reflection operators {t 
the representations are listed. 



, r^}, respectively. Furthermore, the dimensions In of 



APPENDIX C: PROJECTION OPERATORS OF CONTINUOUS SYMMETRIES 



There exists several methods for calculating the projection operators corresponding to a 
continuous symmetry given by a Hermitian operator O; here we think of O being either the 
spin operator S"^ or the pseudospin operator J^. In this Appendix we present an algorithm 
based on successive applications of O. Let |0o) be a given normalized state we want to 
project into 0-symmetry invariant subspaces. When O acts on |0o) a term proportional 
with |(/)o) is generated together with a rest term. The rest term is denoted /i|</>i), and it 
defines a new unit vector perpendicular to |0o) while /i is a prefactor: 

6|0o) = (0o|O|0o)|0o) + /l|0l). (CI) 

If /i = we are done, if not, we continue by applying O to |0i) and expand the result on 
l^o); The rest term is now denoted /2I02), where |02) is a unit vector perpendicular to 
both |0o) and |02) and /2 is a prefactor: 

6|0i) = (0o|O|0i)|0o) + (0i|O|0i)|0i) + /2I02). (C2) 

Since we are working in a finite Hilbert space this process is guarantied to yield a zero rest 
term after M steps, i.e. /a/ = 0. The set >S'q(|0o)) = {|0o), • • • , thus yields the 

smallest 0-invariant subspace containing the starting vector |0o)- The symmetry operator 
O is then diagonalized within 5'q(|0o)) yielding the eigenvalues Uk and eigenvectors \ujk): 

M-l 

d\ujk) = ujk\ujk) with \(^k) = J2 (C3) 

i=0 

The projection V'^'^ of |0o) into the 0-symmetry invariant subspace corresponding to the 
eigenvalue cok is thus simply: 



(C4) 



Based on this equation we find the expressions for the projections in spin space, Eqs. (|TTp- 
(O), and in pseudospin space, Eqs. ([T6|)-(pTD. 



APPENDIX D: SIMULTANEOUS EIGENSTATES OF f AND U 



In this appendix we present the analytical construction of simultaneous eigenstates of T 
and U. The existence of these states explains the remaining degeneracies in the symmetry 
invariant subspaces, and the number of them (found numerically) are listed in Table ^ 
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grouped after the [/-eigenvalue 7 (=0,1,2 for four-electron systems) and the spin S. We 
denote these states T/U states or {ips)- A priori, eigenstates of U is most conveniently 
described in real space whereas eigenstates of T are naturally given in momentum space. 
To require a state to be a T/[/ state imposes severe constraints. Below we find many of 
these states analytically. Since in the following we will be using both real space states and 
momentum space states, we will reserve the letters a, b, c, and d for sites in real space and 
the letters k, q, p and r for momenta. 

We begin from below in Table |T| by studying two-pair states, i.e. 7 = 2. In a naive 
sense, such a state must be a superposition of extremely localized states in real space or 
correspondingly of very out-spread states in momentum space. This involves superpositions 
of many states with different wave vectors k and hence different energies given by Eq. 
It turns out that to obtain an energy eigenstate states only states where k enters together 
with TT — k, where tt = (vr, vr), can be used since cos(/c'') -|- cos(7r — A;^) = 0, and consequently 
these states only exist for even L. We construct the desired state \ipsZo) by superposing 
two-pair states: 

= E e^"-V"-^l«, b; a, 6) = E |k, P; (tt - k), (tt - p)). (Dl) 

ab kp 

From the real space representation it is seen immediately that U\ip'gl^) = 2|-?/'2~q) and that 
'Pf^ \'^'s=o) — 1^5=0)' while the momentum space representation directly yields T\iP1Zq) = 
OlVsIo)- In fact, as can be seen in the rows 7 = 2 of Table |T|, this is the only T/U state with 



7 = 2, e.g. it is easily seen from Eqs. (|TT])-(|T^) that no 7=2 state can have S=l or 2. Note 
how |V'5=o) is independent of the coupling strength U, but that its energy is fZ-dependent 
and of the form E{U) = E{0) + 2U. 



Having explained the 7=2 rows of Table |T| we turn to the 7=1 rows. From Eqs. ([TT|)-(p!3| 



we find that there exists no 7=1 states with S=2. To find the 7=1 states with S'=0,1 we 
construct states iV'kq) which manifestly contains exactly one pair, such that Ul'ipjm) = I'ipkq)'- 

a 

= Y: Ip, k; (tt - p), (q - k)) - 5: |p, r; (tt - p), (q - r)). (D2) 

P ^ pr 



Note how the double sum in Eq. (|D2|) involves all k vectors but that it only depends on q. 



To obtain 5*=! we simply form anti-symmetric combinations of such states: 

l^kq) = IV'kq) - |V'{q-k)q) 

= ^(|p, k; (tt - p), (q - k)) - |p, (q - k); (tt - p), k)). (D3) 
p 

Firstly, we note that I'^/'kq) 7^ if and only if q 7^ 2k. Secondly, since only permutations 
of the vectors k and (q — k) are involved we find that H\tp^^) = E\iIj^^). And thirdly, 

V's^ I'lpkci) = ~l^kq)- Thus we have found ) f/U states lipsZi)- This yields exactly the 
number of states listed in the 7=1/S'=1 row of Table |T|. 

ml) = ICq)- (D4) 
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The states |V's=o) must be sought among hnear combinations of defined as 

l^l^q) = + l^(q-k)q) (D5) 

2 

= H [|P>k; (tt-p), (q-k)) + |p, (q-k); (7r-p), k)]-— ^ |p, r; (7r-p), (q-r)). 

The problem now, however, is that the double sum does not vanish, hence preventing the 
state of being an energy eigenstate. Only upon forming differences iV'kq) — iV'k'q) with 
k' 7^ k can we get rid of it. But the resulting state will only be an energy eigenstate 
if Z]f [cos(/i;'') + cos(g^ — A;'')] equals Z](^[cos(/c"') + cos(g'' — /e"')]. This is easily obtained if 
k' = {q^ — k^, k^), since then we are only permuting the momentum components. 

mZl) = lO - 1^^,)- (D6) 

By accident there can also exist other values of k' which fulfills the requirement, and besides 
combine pairs such as iV'kq) — l^P^'q) '^^^ ^^^o in some cases combine three states, — 
(l^pq) + l^p'q))> or four, (|^kq) + |V'iq))-(|V'pq) + l^p'q)), or eveu more. By a straightforward 
combinatorial search we find the number of states listed in the 7=1/5=0 row of Table |T|, 
and we have thus identified all states in the 7=1 rows, and found them to be independent 
of U but with an energy dependence of the form E{U) = E{0) + U. 

Finally, we turn to the 7=0 rows of Table JTT[ First we note that all S=2 states of the 
systems are found here. This is easily proved by noting that the states with {S, S^) = (2, 0) 
are formed by applying the spin lowering operator (which commutes with H) twice to 
states with (5*, Sz) = (2,2). But the latter states cannot contain any doubly occupied sites 
since that would yield a lower than maximal value of Sz- Clearly, these states as well as 
their energies are independent of U. Then we consider a large class of energy eigenstates 
with = 0, 1 and 7 = 0, which does not require the momentum component vr, and which 
therefore accounts for degeneracies for any value of L. Writing \{k^,k^)) = we 
construct 7=0 states |0) obeying U\(f>) = by symmetrizing one component, say the x 
component, and letting Vg act on the other: 

10) = \k^„ kl- kl k2)sym ® vP\kl kl kl kl), With (D7) 

\h.x 1.x. ix i.x\ — \l.x ux. 1.x l.x\,\ix 1.x ix l.x\,\l.x 1.x. ix l.x\,\i.x 1.x. ix i.x\ 
I "'1 5 "-2 1 "'S 5 "-4 /sym — I "-1 ? "^2 ) "-3 ? "^4 / ' I "'1 5 "-2 i "^4 5 "-3 / ' I "'2 ? "'I ) "-3 5 "-4 / ' I "^2 ? "'I ) "-4 ? "^3 / • 

Since only momentum permutations enter, |0) is clearly an energy eigenstate. The proper 
spin states are found by the standard projections: 

1^1:5') =^ri0)- (D8) 

Simple combinatorics yields 0,30,300 and 1680 S=0 states and 0,90,1050, and 6300 5=1 
states for L = 3,4,5 and 6 respectively. In analogy with the 7=1 case many more T/U 
states can be constructed for L=4 and 6 and 7=0 when the momentum vector vr is taken 
into account. We give one example of a class of such states. For a given momentum vector 
k = (/c^, k^) we define for 6 = x,y, xy the functions nsiX), 

7r^(k) = (r + 7r,F), 7rj,(k) = (r,F + 7r), 7r^y(k) = (/t^ + vr, F + vr), (D9) 
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based on which we introduce two operators H^''^ and H^''^: 

n5''^|ki,k2;k3,k4) = + |ki, k2; ka, k4) + |7r5(ki), 7r5(k2); ks, k4) 

+ |ki,k2;7r5(k3),7r5(k4)) + \ns{ki),ns{k2)-,ns(k3),7is(k4)), ^^^^^ 

n5''^|ki,k2;k3,k4) = + |7r5(ki), k2; 7r5(k3), k4) + |7r5(ki), k2; ka, 7r5(k4)) 

+ |ki,7r5(k2);7r5(k3),k4) + |ki, 7r5(k2); k3, 7r5(k4)). 

Direct inspection shows 17(11^"'^ — n^''^)|ki, k2; k3, k4) = 0, and by enforcing certain con- 
straints on all eight momentum components this state also becomes an energy eigenstate 
with energy E=0, while applying the projector renders the correct spin S=S': 

|VS)=^r(nr-nr)|ki,k2;k3,k4), with k: = TT-kl, n = 1,2,3,4. (Dll) 
Finally, we note that for lattices containing the momentum 7r/2 as is the case for L = 4, 



even more T/U states can be constructed, in accordance with Table III. An example of this 



can be obtained from Eq. ( pil|) . If for example we let 5 = a; then it suffices to enforce the 



constraint = ±7r/2 while allowing any value for the y components. The result is energy 
eigenstates with energy E = J2n'^os{k^). 

We conclude that many of the T/tJ states found numerically have been constructed 
analytically, and in agreement with the numerical findings all these states are independent 
of U, while their energies are of the form E{U) = E{0) + 'yll. The analytic constructions 
reveal that these states are due to a restricted permutation symmetry for the momentum 
components of states in momentum space. 
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